Ctrl+K

Contracting cube

Contracting cube#

In this demo we simulate a unit cube that is fixed at \(x = 0\) and free at \(x = 1\). We use a transversally isotropic material with fiber oriented in the \(x\)-direction.

import dolfin

try:
    from dolfin_adjoint import (
        Constant,
        DirichletBC,
        Expression,
        UnitCubeMesh,
        interpolate,
        Mesh,
    )
except ImportError:
    from dolfin import (
        Constant,
        DirichletBC,
        interpolate,
        Expression,
        UnitCubeMesh,
        Mesh,
    )

import pulse
from fenics_plotly import plot

pulse.iterate.logger.setLevel(10)

# Create mesh
N = 6
mesh = UnitCubeMesh(N, N, N)

# Create subdomains
class Free(dolfin.SubDomain):
    def inside(self, x, on_boundary):
        return x[0] > (1.0 - dolfin.DOLFIN_EPS) and on_boundary


class Fixed(dolfin.SubDomain):
    def inside(self, x, on_boundary):
        return x[0] < dolfin.DOLFIN_EPS and on_boundary

# Create a facet fuction in order to mark the subdomains
ffun = dolfin.MeshFunction("size_t", mesh, 2)
ffun.set_all(0)

# Mark the first subdomain with value 1
fixed = Fixed()
fixed_marker = 1
fixed.mark(ffun, fixed_marker)

# Mark the second subdomain with value 2
free = Free()
free_marker = 2
free.mark(ffun, free_marker)

# Create a cell function (but we are not using it)
cfun = dolfin.MeshFunction("size_t", mesh, 3)
cfun.set_all(0)

# Collect the functions containing the markers
marker_functions = pulse.MarkerFunctions(ffun=ffun, cfun=cfun)

# Create mictrotructure
V_f = dolfin.VectorFunctionSpace(mesh, "CG", 1)

# Fibers
f0 = interpolate(Expression(("1.0", "0.0", "0.0"), degree=1), V_f)
# Sheets
s0 = interpolate(Expression(("0.0", "1.0", "0.0"), degree=1), V_f)
# Fiber-sheet normal
n0 = interpolate(Expression(("0.0", "0.0", "1.0"), degree=1), V_f)

# Collect the mictrotructure
microstructure = pulse.Microstructure(f0=f0, s0=s0, n0=n0)

# Create the geometry
geometry = pulse.Geometry(
    mesh=mesh,
    marker_functions=marker_functions,
    microstructure=microstructure,
)

# Use the default material parameters
material_parameters = pulse.HolzapfelOgden.default_parameters()

# Select model for active contraction
active_model = pulse.ActiveModels.active_strain
# active_model = "active_stress"

# Set the activation
activation = Constant(0.0)

# Create material
material = pulse.HolzapfelOgden(
    active_model=active_model,
    parameters=material_parameters,
    activation=activation,
)

# Make Dirichlet boundary conditions
def dirichlet_bc(W):
    V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
    return DirichletBC(V, Constant((0.0, 0.0, 0.0)), fixed)

# Make Neumann boundary conditions
neumann_bc = pulse.NeumannBC(traction=Constant(0.0), marker=free_marker)

# Collect Boundary Conditions
bcs = pulse.BoundaryConditions(dirichlet=(dirichlet_bc,), neumann=(neumann_bc,))

# Create problem
problem = pulse.MechanicsProblem(geometry, material, bcs)

# Solve problem
pulse.iterate.iterate(problem, activation, 0.1)

2024-03-09 14:07:23,652 [1227] DEBUG    pulse.iterate: Control: [0.0]

2024-03-09 14:07:23,653 [1227] DEBUG    pulse.iterate: Target: [0.1]

2024-03-09 14:07:23,654 [1227] DEBUG    pulse.iterate: Intial number of steps: 5 with step size 0.02

2024-03-09 14:07:23,654 [1227] INFO     pulse.iterate: Iterating to target control (f_21)...

2024-03-09 14:07:23,655 [1227] INFO     pulse.iterate: Current control: f_21 = 0.000

2024-03-09 14:07:23,655 [1227] INFO     pulse.iterate: Target: 0.100

2024-03-09 14:07:23,656 [1227] DEBUG    pulse.iterate: Check target reached: NO

2024-03-09 14:07:23,656 [1227] DEBUG    pulse.iterate: Maximum difference: 1.000e-01

2024-03-09 14:07:23,657 [1227] DEBUG    pulse.iterate: Try new control

2024-03-09 14:07:23,658 [1227] DEBUG    pulse.iterate: Current control: 0.020

2024-03-09 14:07:23,777 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:07:23,778 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:07:23,778 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:07:23,779 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:07:23,779 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:07:23,780 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:07:23,781 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:07:23,781 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:07:23,782 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:07:23,782 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:07:23,783 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:07:24,780 [1227] DEBUG    UFL_LEGACY: Blocks of each mode: 
  99	full

2024-03-09 14:07:24,990 [1227] DEBUG    UFL_LEGACY: Blocks of each mode: 
  18	full

2024-03-09 14:08:29,112 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:08:29,113 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:08:29,114 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:08:29,114 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:08:29,115 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:08:29,115 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:08:29,116 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:08:29,116 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:08:29,117 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:08:29,117 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:08:29,118 [1227] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-03-09 14:08:29,289 [1227] DEBUG    UFL_LEGACY: Blocks of each mode: 
  10	full

2024-03-09 14:08:29,481 [1227] DEBUG    UFL_LEGACY: Blocks of each mode: 
  3	full

2024-03-09 14:08:31,439 [1227] DEBUG    pulse.iterate: Solver did not converge...

2024-03-09 14:08:31,440 [1227] DEBUG    pulse.iterate: 
NOT CONVERGING

2024-03-09 14:08:31,441 [1227] DEBUG    pulse.iterate: Reduce control step

2024-03-09 14:08:31,442 [1227] DEBUG    pulse.iterate: Assign old state

2024-03-09 14:08:31,446 [1227] DEBUG    pulse.iterate: Check target reached: NO

2024-03-09 14:08:31,447 [1227] DEBUG    pulse.iterate: Maximum difference: 1.000e-01

2024-03-09 14:08:31,447 [1227] DEBUG    pulse.iterate: Try new control

2024-03-09 14:08:31,448 [1227] DEBUG    pulse.iterate: Current control: 0.010

*** Warning: PETSc SNES solver diverged in 2 iterations with divergence reason DIVERGED_FNORM_NAN.

2024-03-09 14:08:34,297 [1227] DEBUG    pulse.iterate: 
SUCCESFULL STEP:

2024-03-09 14:08:34,298 [1227] DEBUG    pulse.iterate: Check target reached: NO

2024-03-09 14:08:34,299 [1227] DEBUG    pulse.iterate: Maximum difference: 9.000e-02

2024-03-09 14:08:34,300 [1227] DEBUG    pulse.iterate: Check target reached: NO

2024-03-09 14:08:34,300 [1227] DEBUG    pulse.iterate: Maximum difference: 9.000e-02

2024-03-09 14:08:34,301 [1227] DEBUG    pulse.iterate: Try new control

2024-03-09 14:08:34,301 [1227] DEBUG    pulse.iterate: Current control: 0.020

2024-03-09 14:08:35,747 [1227] DEBUG    pulse.iterate: 
SUCCESFULL STEP:

2024-03-09 14:08:35,748 [1227] DEBUG    pulse.iterate: Check target reached: NO

2024-03-09 14:08:35,748 [1227] DEBUG    pulse.iterate: Maximum difference: 8.000e-02

2024-03-09 14:08:35,749 [1227] DEBUG    pulse.iterate: Adapt step size. New step size: 0.015

2024-03-09 14:08:35,749 [1227] DEBUG    pulse.iterate: Check target reached: NO

2024-03-09 14:08:35,750 [1227] DEBUG    pulse.iterate: Maximum difference: 8.000e-02

2024-03-09 14:08:35,751 [1227] DEBUG    pulse.iterate: Try new control

2024-03-09 14:08:35,751 [1227] DEBUG    pulse.iterate: Current control: 0.035

2024-03-09 14:08:37,187 [1227] DEBUG    pulse.iterate: 
SUCCESFULL STEP:

2024-03-09 14:08:37,188 [1227] DEBUG    pulse.iterate: Check target reached: NO

2024-03-09 14:08:37,188 [1227] DEBUG    pulse.iterate: Maximum difference: 6.500e-02

2024-03-09 14:08:37,188 [1227] DEBUG    pulse.iterate: Adapt step size. New step size: 0.022

2024-03-09 14:08:37,189 [1227] DEBUG    pulse.iterate: Check target reached: NO

2024-03-09 14:08:37,190 [1227] DEBUG    pulse.iterate: Maximum difference: 6.500e-02

2024-03-09 14:08:37,191 [1227] DEBUG    pulse.iterate: Try new control

2024-03-09 14:08:37,192 [1227] DEBUG    pulse.iterate: Current control: 0.058

2024-03-09 14:08:38,928 [1227] DEBUG    pulse.iterate: 
SUCCESFULL STEP:

2024-03-09 14:08:38,928 [1227] DEBUG    pulse.iterate: Check target reached: NO

2024-03-09 14:08:38,929 [1227] DEBUG    pulse.iterate: Maximum difference: 4.250e-02

2024-03-09 14:08:38,929 [1227] DEBUG    pulse.iterate: Adapt step size. New step size: 0.034

2024-03-09 14:08:38,931 [1227] DEBUG    pulse.iterate: Check target reached: NO

2024-03-09 14:08:38,931 [1227] DEBUG    pulse.iterate: Maximum difference: 4.250e-02

2024-03-09 14:08:38,932 [1227] DEBUG    pulse.iterate: Try new control

2024-03-09 14:08:38,932 [1227] DEBUG    pulse.iterate: Current control: 0.091

2024-03-09 14:08:40,666 [1227] DEBUG    pulse.iterate: 
SUCCESFULL STEP:

2024-03-09 14:08:40,667 [1227] DEBUG    pulse.iterate: Check target reached: NO

2024-03-09 14:08:40,667 [1227] DEBUG    pulse.iterate: Maximum difference: 8.750e-03

2024-03-09 14:08:40,668 [1227] DEBUG    pulse.iterate: Adapt step size. New step size: 0.051

2024-03-09 14:08:40,669 [1227] DEBUG    pulse.iterate: Check target reached: NO

2024-03-09 14:08:40,669 [1227] DEBUG    pulse.iterate: Maximum difference: 8.750e-03

2024-03-09 14:08:40,670 [1227] DEBUG    pulse.iterate: Change step size for final iteration

2024-03-09 14:08:40,670 [1227] DEBUG    pulse.iterate: Try new control

2024-03-09 14:08:40,671 [1227] DEBUG    pulse.iterate: Current control: 0.100

2024-03-09 14:08:41,554 [1227] DEBUG    pulse.iterate: 
SUCCESFULL STEP:

2024-03-09 14:08:41,555 [1227] DEBUG    pulse.iterate: Check target reached: YES!

2024-03-09 14:08:41,556 [1227] DEBUG    pulse.iterate: Check target reached: YES!

([Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 47),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 175),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 224),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 273),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 330),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 387),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 420)],
 [Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 45),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 173),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 222),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 271),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 328),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 385),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 418)])

# Get displacement and hydrostatic pressure
u, p = problem.state.split(deepcopy=True)

V = dolfin.VectorFunctionSpace(mesh, "CG", 1)
u_int = interpolate(u, V)
new_mesh = Mesh(mesh)
dolfin.ALE.move(new_mesh, u_int)

fig = plot(mesh, show=False)
fig.add_plot(plot(new_mesh, color="red", show=False))
fig.show()